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Development of a Linear Stirling System Model 
With Varying Heat Inputs 

Timothy F. Regan and Edward J. Lewandowski 
Sest, Inc. 

Middleburg Heights, Ohio 44130 


Abstract 

The linear model of the Stirling system developed by NASA Glenn Research Center (GRC) has been extended to 
include a user-specified heat input. Previously developed linear models were limited to the Stirling convertor and 
electrical load. They represented the thermodynamic cycle with pressure factors that remained constant. The 
numerical values of the pressure factors were generated by linearizing GRC’s nonlinear System Dynamic Model 
(SDM) of the convertor at a chosen operating point. The pressure factors were fixed for that operating point, thus, 
the model lost accuracy if a transition to a different operating point were simulated. Although the previous linear 
model was used in developing controllers that manipulated current, voltage, and piston position, it could not be used 
in the development of control algorithms that regulated hot-end temperature. This basic model was extended to 
include the thermal dynamics associated with a hot-end temperature that varies over time in response to external 
changes as well as to changes in the Stirling cycle. The linear model described herein includes not only dynamics of 
the piston, displacer, gas, and electrical circuit, but also the transient effects of the heater head thermal inertia. The 
linear version algebraically couples two separate linear dynamic models, one model of the Stirling convertor and one 
model of the thermal system, through the pressure factors. The thermal system model includes heat flow of heat 
transfer fluid, insulation loss, and temperature drops from the heat source to the Stirling convertor expansion space. 
The linear model was compared to a nonlinear model, and performance was very similar. The resulting linear model 
can be implemented in a variety of computing environments, and is suitable for analysis with classical and state 
space controls analysis techniques. 


Nomenclature 


Ad 

Displacer area (m 2 ) 

Qin 

Heat delivered by the heat source (W) 

A P 

Piston area (m 2 ) 

Qout 

Heat rejected (W) 

^ rod 

Displacer rod area (m 2 ) 

Pall 

Alternator resistance (Q) 

C h 

Heater head thermal capacitance (J/K) 

R c 

Compression space resistance (W/m-K) 

C, 

Tuning capacitance (pF) 

Re 

Expansion space resistance (W/m-K) 

D d 

Displacer damping (N-s/m) 

Rgas 

Gas constant (J/kg-K) 

D P 

Piston damping (N-s/m) 

R h 

Heater head resistance (W/m-K) 

f 

Operating frequency (Hz) 

Rhh 

Heater head conduction loss resistance (W/m- 

I alt 

Alternator current (A) 


K) 

R bounce 

Bounce space spring rate (N/m) 

Rins 

Insulation resistance (W/m-K) 

K d 

Displacer spring rate (N/m) 

t 

time (sec) 

R magnet 

Alternator magnet space spring rate (N/m) 

T a 

Ambient temperature (K) 

K P 

Piston spring rate (N/m) 

T a lt 

Alternator temperature (K) 

L a /t 

Alternator inductance (H) 

T 

X C 

Compression space temperature (K) 

M d 

Effective displacer mass (kg) 

T e 

Expansion space temperature (K) 

M p 

Effective piston mass (kg) 

T h 

Hot heat exchanger temperature (K) 

MW 

Gas molecular weight 

T k 

Cold heat exchanger temperature (K) 

N 

Number of turns on the alternator winding 

T 

x r 

Regenerator space temperature (K) 

P 

Stirling cycle dynamic pressure (Pa) 

T 

x source 

Heat source temperature (K) 

Pc 

Compression space pressure (Pa) 

v c 

Compression space volume (m 3 ) 

Pd 

Displacer pressure factor (Pa/m) 

v co 

Mean compression space volume (m 3 ) 

Pc 

Expansion space pressure (Pa) 

v e 

Expansion space volume (m 3 ) 

p 

1 exp 

Expansion space PV power (W) 

Veo 

Mean expansion space volume (m 3 ) 

P P 

Piston pressure factor (Pa/m) 

v h 

Hot heat exchanger volume (m 3 ) 

Qe 

Heat into the Stirling cycle (W) 

v k 

Cold heat exchanger volume (m 3 ) 
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V T 

Regenerator volume (m 3 ) 

mag 

Alternator magnetic efficiency 

x d 

Displacer position (m) 

a> 

Flux (Wb) 

x d 

Displacer position amplitude (m) 

TC 

Displacer phase angle (deg) 

Xp 

Piston position (m) 

% 

Pressure phase angle (deg) 


Piston position amplitude (m) 

CO 

Operating frequency (rad/sec) 


AP Pressure drop across regenerator, heater and 
cooler (Pa) 


I. Introduction 

NASA’s Glenn Research Center is using its System Dynamic Model (SDM) to develop linear models of free- 
piston Stirling convertors for use in real-time simulations of Stirling convertors. The outputs of such a simulation 
could be used in place of actual convertor hardware in tests of controllers and other convertor support hardware. The 
SDM is a non-linear model that simulates Stirling convertor system dynamics for systems of arbitrary complexity. 
SDM has been used to simulate a dual-convertor system with controllers based on pulse-width modulated (pwm) 
power electronics, however the maximum time step allowed for a 50 kHz pwm controller simulation is on the order 
of nano-seconds. Consequently, no more than a few cycles of such a system can be simulated on SDM running on a 
typical desk-top computer. On the other end of the spectrum, when the slow variation of system temperatures 
requires study, the simulation times can become so long as to make the simulation approach impractical. To address 
these extremes, SDM was used to linearize the model at an operating point and the linearized coefficients were used 
to build a linear model. Linear models run faster than SDM and they can run in a variety of computing hardware 
environments. Linear models of Stirling convertor systems, including a dual-opposed system and a single convertor 
system with a balancer were analyzed in detail as part of NASA GRC’s Stirling convertor modeling effort 1 . The 
models reported in the literature have included only the mechanical and mounting dynamics. The pressure wave was 
modeled by constants — called pressure factors — that were calculated by SDM to represent the pressure wave at the 
point at which the model was linearized. The models described in Ref. 1 were not equipped to track the relationship 
between output power, hot-end temperature and piston amplitude. Convertor conditions which cause an increase or 
decrease of hot-end temperature such as a change in piston amplitude, a change in the input heat or the rejection 
temperature were not tracked successfully by those linear models. The models lost some accuracy whenever the 
operating point moved away from the operating point at which the model was linearized. That is, the more the piston 
amplitude, hot-end temperature or output power changed during a simulation using the models described in Ref. 1, 
the less accurate the results became. Such a model was not useful in controlling against disturbances in thermal 
environment or user load since disturbances in these areas invalidate the model. 

The linear models developed so far have followed the state-variable formulation in which the system matrix — 
the A-matrix — was constrained to have constant coefficients. If the system changes with operating conditions, then 
the coefficients in the model cannot remain constant. In the present work, a method is described for re-computing 
the hot-end temperature based on the calculated expansion space PV power and re-computing certain entries in the 
A-matrix as the system temperatures change. 

The validation data for the models described herein were published by NASA in 1999 from operation of the 
Component Test Power Convertor (CTPC). The CTPC effort published design details as well as test data. The final 
report 2,3 has been used for validation of Stirling simulation codes and design codes. A model of the CTPC was 
constructed in SDM and used to produce linearized A-matrix coefficients as in the models described in Ref. 1. 
Section II begins with a description of the CTPC and its SDM model. Section 1I-A describes the linear electro- 
mechanical model. Section II-B describes the thermal system model. Section II-C describes how the models are 
coupled through re-computation of the pressure factors. Section III summarizes the results. 

II. CTPC Dynamic Model 

A high-level schematic of the SDM model of the CTPC is shown in figure 1 . Starting at the bottom of the figure, 
the heat input to the Stirling cycle was represented by the heat input Q in , shown with two flow paths emanating from 
it. One flow path represents the heat lost through the insulation to ambient and has the thermal resistance R ins . The 
other flow path represents the heat conducted to the heater. There is a temperature drop from the heat source to the 
heater. The temperature drop was modeled by the resistance R/,. The thermal time constant of the heat input system 
was modeled by the thermal capacitance C h . This parameter was determined based on the mass of the CTPC heater 
and its heat capacity. 
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Not all of the heat entering the heater in the Stirling convertor 
goes into the Stirling cycle. Some is conducted to the cold end 
through paths including the heater head wall, inner cylinder, 
helium gas, regenerator matrix and the displacer. This loss was 
modeled by the heat flow through the thermal resistance R/,/,. 

The Stirling cycle thermodynamics were modeled based on 
the Schmidt model 4 , which assumes an isothermal Stirling cycle, 
but in which pumping losses through heat exchangers and the 
regenerator were considered. The thermodynamic portion of the 
model determines the internal expansion space and compression 
space gas pressures P e and P c . These pressures generate the 
driving forces on the displacer mass M d and the piston mass M p . 
In the model, gas spring and damping forces on the displacer were 
represented by K d and D d . Displacer damping is a non-linear 
function of fluid flow through the porous regenerator material and 
through the heat exchanger tubes. The piston is subjected to 
conservative forces from the bounce space gas spring and the 
magnet spring. In the model, these were represented by K bounce 
and K magnel . The electrical load circuit also contributes spring 
force if the reactance of the alternator winding is not completely 
balanced by the capacitive reactance of the tuning capacitor. . 

The alternator produces a damping force on the piston based 
on the alternator current. Alternator electrical dynamics were 
determined by the alternator resistance, inductance, motor 
constant, and the tuning capacitor. 

The Stirling convertor was assumed to be rigidly connected 
to the ground. Casing motion, dual-opposed dynamics, and the 
effect of a dynamic balancer could be added to this model if 
desired. 



Figure 1. — Schematic of the CTPC 
dynamic model. 


A, Linear Dynamic Stirling Model 

The ultimate design of a model is determined by a number of factors, including fidelity required, simulation 
time, and model platform. Linearized parameters were used to create a model in Simplorer (Ansoft Corporation) and 
other dynamic modeling platforms. The equations for the electro-mechanical portion of the model were developed in 
Reference 1 and are reprinted as eqs. (1) through (4) below:. 




dt 


= x d 


dx d 

dt 


f „ , dP 

K d + A 


V 


8x d ) M d 


x d + 


r A SAP'] 1 . , 31> 1 

Afj _ | A r _ _ _ x n + 


V 


dx d JM d 
dx„ 


8x„ M d 


SAP 
8x p ) 




dt 


- = x. 


( 2 ) 


( 3 ) 


dx p 8P 1 

— = -A d x d - 

P 8x d M p 


dt 


K„ + A, 


8P 


^ p 8x 


p 


1 °p . M 

x„ — x„ +N 

M p P M p p 


d<P 1 1_ 

dx p q mag Mp 


I alt 


( 4 ) 


Equations (1) through (4) represent the dynamic equations for the displacer and the piston. In eq. (2), the force 
due to the pressure wave’s action on the displacer area A d has been accounted for in two parts. This is because the 
pressure is expressed as a linear combination of the piston and displacer positions. To reconstruct the pressure wave, 
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the contribution of both components must be considered. The same discussion applies to the pressure force on the 
piston in eq. (4). In state-variable format, 


X = A-X + Bu 


(5a) 


with the state vector defined as 


X = [x d 


Xd 



( 5 ) 


the ^-matrix entries are seen to be the coefficients of the state variables in eqs. (1) through (4), above. In eq. (4), the 
right-most term expresses the piston damping due to alternator current. In a state-variable formulation, I a i, is the 
input, u, and the coefficient of I ait is the 5-matrix entry that contributes to eq. (4). Thus, 


B = 


0 0 0 N 


d <. P 1 

dx p T"| ma g 



(6) 


B. Thermal System Model 

Figure 2 shows the thermal system modeled as a series of 
thermal resistances with a heat capacitance, heat input, and 
temperature sinks. It is analogous to figure 1 with the 
electromechanical components removed. The schematic 
shows the system as it is configured in the Simplorer 
modeling platform 

Thermal energy from the heater Q in is shown to be 
divided between two flow paths: 1) heat lost through the 
insulation to ambient was determined from R ins and 2) heat 
transferred to the heater head was determined from R/,. The 
resulting heat flow to the heater head causes an increased 
convertor hot-end temperature T h . 

Thermal energy from T h is shown to be divided among 
three flow paths: 1) heat lost via conduction through the 
convertor was determined from Ru„ 2) heat transferred to the 
gas in the expansion space is determined from R e , and 3) heat Figure 2.— The thermal system was modeled as a 

transferred to the heater head material was determined from circuit in SDM. 
the heater head thermal capacitance C/,. Thermal capacitance 

only becomes apparent during transient analysis; at steady-state operation the heater head material temperature 
equals 7/„ and heat flow is divided between the conduction loss and heating of the gas expansion space, with the 
bulk of the heat being transferred to the latter. 

Heat is rejected from the convertor to the cold-sink temperature T k . Most of the rejected heat flows from the 
working fluid. The remainder of the heat rejected is from the conduction losses. The circuit model shown in figure 2 
displays the operation of the thermal system. Circuit solvers such as the one included with Simplorer offer a 
convenient way to calculate the temperatures and heat flows throughout the thermal system. Simplorer will not 
simulate the solution in real time however. To implement a Stirling convertor simulator in real time, the circuit 
solution must be implemented in more direct mathematical language. The method presented here follows a state 
variable formulation. There is one energy storage element in the thermal system: the heater head thermal mass, Q, 
therefore, there is only one state variable in the formulation. 

1. State Variable Formulation of the Thermal System 

In the circuit diagram of figure 3, as in the circuit of figure 2, there is one energy storage element, the heater 
head thermal mass, C/,. The single state variable is the temperature across it, T h . There are four inputs — the heat in, 
the expansion space PV power, the ambient temperature, and the rejection temperature. The inputs are Q,„, Q pv , T a , 
and T k , respectively. Heat source Q m models a fixed heat source such as a radioisotope heat source. Q pv is a 


Qin, 

Tsource Fh Th 
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Figure 3. — A simplified circuit was used to derive the state-variable equation 
used in the linear model. 

continuous calculation of the expansion space work over the cycle. Observe that when Q in is greater than Q pv , the 
excess heat flows into the capacitor, C h . The temperature across C/, will rise as a result. The heat through a thermal 
mass (capacitance) is related to the temperature across it as shown in eq. (7). 

Qc=Qin-Qpv=C h d ^. (7) 

If the heat flow into the capacitor increases, then Q in will be greater than Q pn and the temperature across the 
capacitor will increase. If the heat flow into the capacitor decreases, then Q, n is less than Q pv , and the temperature 
across the capacitor will decrease. 

The circuit of figure 3 is more simple than the circuit of figure 2 in that the thermal resistor, R h has been 
eliminated. This was done to simplify the math used to calculate the expansion space temperature. The process of 
reducing the thermal circuit of figure 3 to state-variable form was motivated by the observation that the heat 
supplied to the convertor by (Q,„ is the sum of the insulation loss, the heater head conduction, the expansion space 
PV power, and the heat into the thermal capacitance 


dT 

Qin ~ Qins Qhh ^ Qpv Qc ~ Qins "t" Qhh ^ Qpv ^ ii I 


dt 


( 8 ) 


Additionally, the insulation loss and the heater head conduction can be expressed in terms of the inputs and the 
state variable as shown in eqs. (9) and (10). 


Qins 



(9) 


Qhh = 


Th ~Tk 

R-hh 


(10) 


The expansion and compression space temperatures, T e and T c can be expressed in terms of input, Q pv and the 
state variable, 7)„ thus: 


T e =T h — R e ' Qpv and T c = T h + R c -Q pv 

By substituting eqs. (9) and (10) into eq. (8) and simplifying, the following expression was obtained: 

dTh _ ~Th Qhh 1 + R ins 1 ) | Qin QpV T k T A 

dt C h C h C h C h • Rhh C h • R j ns 

Equation (12) is in the basic form of the state equation. 


( 11 ) 


( 12 ) 
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X = AX + Bu 


(13) 


where A is the single coefficient 


A — ci 11 


- (R hh 1 + R ins 1 ) 
C h 


(14) 


and B is the vector [b] b 2 bj h.,] 1 where 



(15a) 


bi 



(15b) 


£3 


1 


Ch ■ Rh 


(15c) 


b 4 = 


1 

C h ' Rhh 


(15d) 


Thus, there is a B-vector coefficient for each of the four inputs, 
Qin, Qpv, Ta and T k . When the product Bu is evaluated, one number is 
produced. It has dimensions of Kelvins per second. An 
implementation of the state equation in block diagram form is shown 
in figure 4. 

The thermal model determines the hot-end temperature (T h ) for 
the Stirling cycle, but does not calculate cycle performance or 
convert heat energy into mechanical energy. This is accomplished by 
coupling the linear Stirling model to the thermal system model. The 
next section will discuss how the model coupling was performed. 


C. Coupling Stirling Linear Model and Thermal System 

Model 

To complete the linear model of the CTPC, the Stirling dynamic 
model has to be coupled to the thermal system model. One approach would be to expand the Stirling cycle model of 
eqs. (1) to (4) to include the thermal system. The concern with this approach is that the response time of the Stirling 
cycle is on the order of milliseconds, while the response time of the thermal system are on the order of seconds or 
minutes. By keeping the electro-mechanical system separate from the thermal system, the likelihood of formulating 
an ill-conditioned matrix is reduced. 

An alternative approach is to algebraically couple two separate linear dynamic models, one model of the 
Stirling cycle and one model of the thermal system, through the pressure factors. The pressure factors are the partial 
derivatives of the pressure wave P with respect to displacer position x d and piston position x p . The equations to 
couple the systems can be derived from the Stirling cycle equations as derived by Urieli and Berchowitz 4 . 



b1*Qin+b2*Qpv+b3*Ta+b4*Tk 

Figure 4. — Thermal model block diagram 
used to solve the state equation. 


P = My., ■ R 


gas 


YIl + Yl. 
T,. 






\-l 


(16) 
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In eq. (16) the volumes V e and V c are functions of x d and x p . The other variables are assumed to be fixed over a given 
cycle. 


V e =V eo -A d x d (17) 

Vc — V c0 — ApXp + (A d — A rod )x d (18) 

The temperatures T c and T e can be expressed in terms of the heat into the heater head Q k based on the thermal circuit 
shown in figure 2, where the thermal resistance R e is used to calculate the temperature drop from T h to T e , and R, is 
used to calculate the temperature drop from T c to T k . 

T e - T h ~ Qh R e (19) 

T c = Q h R c + T k (20) 


The regenerator temperature can be expressed as the log mean temperature between T h and T k : 

j T h ~ Tk 

r HT h /T k ) 


(21) 


Substituting eqs. (17) to (20) into eq. (16) and differentiating yields the following expressions for the pressure 
factors 


dP 

8x n 


A p M w R gas 
{Q e R c + T k ) 


Vco dp x p + (^d A-od ) x d + Vk_ + V r In (Th / T k ) + Yh_ + Veo ~ Al x d 
Q e R c + T k T k T h - T k T h T h - Q e R e 


1-2 


(22) 


8P --M R 

A ro d ) 

A d 

Vco A p Xp+(A d A wd )x d | V k | V r \n(T h !T k ) ^ V k | V eo -A d x d 

— Ivl w^gas 

dx d 

_ QeRc + Tk 

{Th — Q e R e ) _ 

Q e R c + T k T k T h - T k T h T h - Q e R e 


1-2 


(23) 


The pressure factors are now functions of T h , T k , Q e , which can be determined based on the dynamics of the thermal 
system. The pressure factors are also functions of x p and x d . In theory, the pressure factors should be the same 
regardless of the point in the cycle at which x p and x d were measured, but analysis of simulation data has shown that 
the calculated values of the pressure factors do vary over the cycle. It has been observed that the calculated values of 
the pressure factors are closest to the actual value as determined from analysis of the pressure wave when x p = x d and 
are both negative. Further analysis is needed to understand the reason for this and to understand the limitations. 
Given displacer and piston amplitudes X d and X p and a displacer phase angle (p,/, x p and x d are equal and negative 
when 


x p = -X p sin 


tan 


-l 


X d sinf/ 

X - X d cos(i> rf 


/ 

Xd = -Xd sin 

v 


tan -1 


x d sinjy 
Xp -X d coscl),/ 


V 

+ ( t ) rf 

/ 


(24) 


(25) 


These values for x p and x d should be used in eqs. (22) and (23). 

To algebraically couple the thermal system with the Stirling model, it was also necessary to calculate the 
expansion space PV power. The expansion space PV power is approximately equal to the heat flow into the Stirling 
cycle 
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For the expansion space, 


Q h * ffrdV 


(26) 


^eo '^t/ X d 


Since 


x d (t) = X d sin(( 0 l + <j> rf ) 


(27) 


(28) 


Xp(t) = X p sin(cot) 


(29) 


then 


— - = -A d X d o)cos(o)t + <|> rf ) 
dt 

Also, the pressure in the expansion space is a function of the pressure factors 

-P(t) = P d X d ( t ) + PpXp (t) 

Substituting eqs. (30) and (31) into eq. (26) and integrating gives the PV power: 

Q P v = KA df X d XpP p sin <j)j 


(30) 


(31) 


(32) 


Determining the PV power is straightforward by use of eq. (32) if the phase angle and amplitudes are known. 
SDM makes use of Simplorer’s state diagrams to calculate them from the piston and displacer position signals. The 
same approach can be taken with Matlab/Simulink (The MathWorks, Inc.) using the Stateflow (The MathWorks, 
Inc.) add-on package. Implementations that make use of the Very-high-speed integrated circuit Hardware Definition 
Language (VHDL) are able to use the state machine definition features of VHDL to determine amplitudes and 
angles. The amplitude of a sinusoidal signal can be determined by observing when the derivative of the signal 
crosses zero. Since the derivative of a sinusoid is a sinusoid lagging by 90°, the zero-crossing of the derivative 
coincides with the maximum or minimum of the original signal. Detection of zero-crossings can be detected by state 
machines programmed to change states when the control signal changes signs. 

An even more direct solution was used in the Simplorer implementation of the linear model that did not require a 
state machine. The PV power was calculated by the product of frequency and the output of a power measurement 
block. The two inputs to the power measurement block were expansion space volume and the pressure wave instead 
of voltage and current. 

If temperature changes in the hot end, then the working fluid redistributes itself within the working space and the 
non-working spaces of the convertor. The non-working spaces comprise the bounce space and the interior of the 
displacer. The mass of fluid in the working space is computed from the total charge, M tot , as follows: 


+ V r + y k + V e< 


V„ 


M W =M 


tot ' / 


v Th Tr T k T e T, 


c j 


V h + K + v L+ v i , 


l T h 


Pgo ^ bounce 

T T 

1 c 1 bounce y 


(33) 


The pressure factors as calculated by eqs. (22) and (23) were used to re-compute the entries for a 2 1 , a 2 3 , and a 4 1 
in eqs. (34) through (36). 
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f 


a 2\ ~ ' 


K r , +A 


dP 1 


rod 


a 23 ~ ~ 


dP 


dx d ) M d 
\ 


*rod 


dx 


1 


P 


M n 


a 4\ - _ 


dP ) 1 


p dx d ) M, 


a 43 - ■ 


K „ + A n ■ 


dP 


dx 


p 


M r 


(34) 

(35) 

(36) 

(37) 


The only other coefficients to be calculated are 
those responsible for the displacer damping. Their 
values depend on the pressure drop between the 
expansion and compression spaces. The pressure drop 
was modeled as a sinusoidal pressure expressed as a 
linear combination of the piston velocity and the 
displacer velocity. Values for these pressure drop 
factors were determined by performing an SDM model 
simulation through a temperature range and fitting the 
model’s recorded output to a linear curve with hot-end 
temperature as a parameter. 


xl 



Figure 5. — Simplorer implementation of 4-state, 2-input 
state equation uses equation block and four integrators. 
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Finally, a 2 2 and a 2 4 were calculated: 
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D. Implementation 

The linear models and their coupling equations were implemented in Simplorer and in Matlab/Simulink. These 
two high-level modeling platforms both offer convenient means to compute the integral of the state equations and to 
evaluate the output power and the expansion space PV power. In Simplorer, the integration of a state equation is 
accomplished with the block diagram of figure 5. Simulink offers a ready-made block for state equations which 
requires only the specification of a 4 x 4 system matrix (A-matrix) and a 4 x 2 input matrix (B-matrix) to implement 
the same system shown in figure 5. Implementation of the linear model also required the equations for re-computing 
the six A-matrix entries, and evaluating the expansion space PV power. 
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Figure 6. — Notional convertor simulator is based on the linear model running in the FPGA and 
producing an EMF signal in real time. 

There are a number of software tools available for translating Matlab/Simulink models into code that can be 
executed on a microprocessor or a field-programmable gate array 6 (FPGA). The advantage of running the model on 
such hardware rather than on a desktop computer is that the model can simulate convertor performance in real time. 
With real-time simulation, the actual thermal time constants can be used in the model and more realistic response 
obtained. Simplorer offers a feature whereby the model shown in figure 4 can be translated into VFIDL code. The 
feature can be used to generate a first cut of the VHDL code used to program the FPGA to execute the linear model. 
The code generated is not directly synthesizable — the user must first substitute the non-synthesizable variable 
declarations with appropriate fixed-point declarations. 

By implementing the linear model on a readily-available demo board for a microprocessor or FPGA, and connecting 
the EMF output to a power amplifier, a convertor simulator such as that shown in figure 6 can be constructed. 

The State variable representing the piston velocity is multiplied by the alternator constant and the result is used 
as the main output of the model. That output is equal to the alternator EMF. It is amplified in the power amplifier 
and the output is wired in series with an inductance and a resistance. The values used are equal to the inductance and 
resistance of the alternator winding. The power measured by the watt meter is the actual power output of the 
modeled convertor. In the case of the CTPC, it is 12.5 kW, so the amplifier in the notional schematic of figure 6 is 
not one that is readily available. The voltage on the load end of R a /, is the same as the voltage at the terminals of the 
convertor. Controller hardware and software can be connected and exercised. Various parameters can be changed in 
such a system, perhaps even while the system is running. The changes can be introduced during a controller test to 
verify that the controller performs as expected when the convertors endure radiator temperature changes as well as 
the expected radioisotope decay. The real-time simulator based on the linear model could also find use as a system 
for demonstrating design margins in convertor matching for dual-opposed pairs. For example, two simulators and 
two controllers could be used to compare system performance with spring rates matched to within 0.1, 1, and 10%. 

III. Summary 

The linear model described extends the four-state electro-mechanical linear model previously developed from 
outputs of SDM. The extended model includes a user-specified heat input and a hot-end temperature that varies with 
operating conditions. 

Seven mathematical constructs were required in order to implement the model 

1. A four-state, single-input, single-output state space model of the piston and displacer motion and 
damping as shown in eqs. (1) through (4). 

2. A single-state, four-input, single-output state space model of the thermal processes in the hot end as 
shown in eqs. (12) and (13). 

3. Recalculation of the mass of the working fluid in the working spaces as shown in eq. (33) whenever a 
change in temperature is detected. 

4. Recalculation of the pressure factor of eq. (20) and associated A-matrix entries, a 2 3 (eq. (33)) and fl 43 
(eq. (35)) whenever a change in either temperature, amplitude or phase angle is detected. 

5. Recalculation of the pressure factor of eq. (21) and associated A-matrix entries, a 2 1 (eq. (32)) and « 4 | 
(eq. (34)) whenever a change in either temperature, amplitude or phase angle is detected. 

6. Recalculation of the pressure drop factor for displacer damping with respect to displacer velocity from 
eq. (36) and the associated A-matrix entry a 22 from eq. (38). 

7. Recalculation of the pressure drop factor for displacer damping with respect to piston velocity from 
eq. (37) and the associated A-matrix entry a 2 4 from eq. (39). 
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A possible use of the model is in Stirling convertor simulators. The simulator will generate realistic convertor 
electrical terminal voltages and currents so that convertor systems including controllers can be exercised under a 
variety of conditions. 
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